Pseudo-High-Order Symplectic Integrators 

J. E. Chambers 

Armagh Observatory, Armagh, Northern Ireland, BT61 9DG, United Kingdom 

jec @star. arm. ac.uk 

M. A. Murison 
Q^ , Astronomical Applications Dept., U.S. Naval Observatory, 3450 Mass. Ave., NW, Washington 



o 

O 



% 



D.C. 20392-5420 
murison@aa. usno.navy. mil 

ABSTRACT 



Symplectic N-body integrators are widely used to study problems in celestial me- 

chanics. The most popular algorithms are of 2nd and 4th order, requiring 2 and 6 

^r^ • substeps per timestep, respectively. The number of substeps increases rapidly with 

^«0 . order in timestep, rendering higher-order methods impractical. However, symplectic 

f^ . integrators are often applied to systems in which perturbations between bodies are a 

small factor e of the force due to a dominant central mass. In this case, it is possible to 

Qv^ [ create optimized symplectic algorithms that require fewer substeps per timestep. This 

(— I I is achieved by only considering error terms of order e, and neglecting those of order 

Qh' e^, e^ etc. Here we devise symplectic algorithms with 4 and 6 substeps per step which 

O ' effectively behave as 4th and 6th-order integrators when e is small. These algorithms 

-4— > ' are more efficient than the usual 2nd and 4th-order methods when applied to planetary 

C^ ■ systems. 



Subject headings: celestial mechanics, stellar dynamics — methods: n-body simulations- 
methods: numerical 



1. Introduction 

Symplectic integrators are widely used to study problems in celestial mechanics. These inte- 
grators have two advantages over most other algorithms. First, they exhibit no long-term build up 
in energy error. Second, the motion of each object about the central body can be "built in" , so that 
the choice of step size, r, is determined by the perturbations between bodies, whose magnitude is 
a factor e smaller than the forces due to the central body ((Wisdom and Holman 1991)). 

The most popular algorithm is the second-order symplectic integrator. The error at each step 
is proportional to er^, so that the likely error for an integration as a whole is ~ er^. The second- 
order method is easy to implement, consisting of only two substeps, including one force evaluation, 
per time step. It is also very fast for integrations requiring moderate accuracy. 



For more accurate integrations, it is better to use the fourth-order method ((Forest and Ruth 1990)). 
Here, the error at each step is proportional to er^, although each step is computationally more ex- 
pensive since it consists of 6 substeps. Yoshida (1990) has developed 6th and 8th-order symplectic 
integrators. However, these do not appear to be competitive in most situations, due to the large 
number of substeps required. 

Here we show how to construct what are effectively high-order (4th, 6th etc.) symplectic 
integrators that require fewer substeps per time step than those in current use. The trick is to take 
into account the dependence of each error term on e when choosing the coefficients for each substep. 
The algorithms are designed by eliminating error terms proportional to e up to the desired order 
of the timestep. Error terms proportional to e^,e^ etc., in low orders of the timestep, still exist. 
However, in many situations these terms are negligible, and the integrators behave as if they are 
of higher order than the leading error term in r suggests. 

Section 2 gives a quick review of how symplectic integrators are traditionally constructed using 
Lie algebra. In Section 3, we show how to build more efficient symplectic algorithms using fewer 
substeps. Section 4 contains results of test integrations that compare the new algorithms with 
traditional symplectic integrators. 



2. Symplectic Integrators 

Symplectic integrators for the N-body problem can be constructed starting from Hamilton's 
equations of motion: 

dxi dH 

dt dpi 

dpi _ _dH^ , . 

dt ~ dxi ^ ' 

where Xi and pi are the coordinates and momenta of each body respectively, and H is the Hamil- 
tonian of the system. 

Using these equations, the rate of change of any dynamical quantity g(x, p, t) can be expressed 
as 

dq y--^ / dq dH dq dH\ _ r rj^ _ „ 
dt ^ \ dxi dpi dpi dxi J 

where {,} are Poisson brackets, and F is a differential operator. 

The formal solution of equation (2) is 



(t) = e-^ q(^t-T)=[l + TF+ -^- + ...]q{t-T) 



where F'^q = F{Fq) etc. 



Now suppose that we are able to split the Hamiltonian into two pieces, Ha and Hb, so that 
each part of the problem can be solved relatively easily in the absence of the other. The solution 
for q becomes 

q{t) = e^^^+^U{t - t) (3) 

where A and B are differential operators related to Ha and Hb respectively, in the same way that 
F is related to H. 

The Baker-Campbell-Hausdorff (BCH) formula states that, for any noncommutative operators 
A and B, 

exp{A) • exp(i?) = exp(C) 

where C is a series consisting of nested commutators, 

C = A + B + ^[A,B] + ^[A,A,B] + ^[B,B,A] + ... 

where the commutator [A, B] = AB — BA 7^ in general (see, for example, Yoshida 1990 or Forest 
and Ruth 1990). Here, we have used the nested commutator notation [A, i3,C] = [A, [B,C]], etc.. 

Hence, if we evolve q under the two parts of the Hamiltonian separately, one after the other, 
we have 



exp(rA) • exp(ri?)g'(t — r) = exp 



r2 



rF + -[A,B] + ... 



q{t - t) (4) 



This is identical to the righthand side of equation (3) to 0{t), and so equation (4) represents a 
first-order integrator. Each step of the integrator consists of 2 substeps, with the whole step giving 
an error of 0{t'^). Alternatively, we can say that the integrator exactly solves a problem whose 
Hamiltonian is given by 



Hinteg = H + -{Hb,Ha} + ©(t^^ 



T 
2' 

(see, for example, (Saha and Tremaine 1992)). Provided that r is small, and {Hb,Ha} remains 
bounded, the energy of the integrated system will always be near to that of the real system. 



Other integrators can be found by combining exponential operators in such a way that they 
are equivalent to equation (3) up to a given order in r. For example, we have the second-order 
symplectic integrator 

S2A = exp (-A) ■ exp(r5) • exp (-A 



exp 



r3_ „ ., r3 



rF + -[B,B,A]--[A,A,B] + ... 



When many integration steps are performed one after another, the exp(r^/2) terms at the end of 
one step and the start of another can be combined. Hence, the second-order integrator also consists 
of only 2 substeps, except at the beginning and the end of an integration. 



Another second-order integrator is 

'T 



S2B 



exp 1 -Bj ■ exp(r^) • exp l^-B 



exp 



3 3 

rF+^[A,A,B]-^[B,B,A] + ... 



(5) 



The distinction between S2A and S2B (which at first sight appear to be the same) will become 
apparent in the next section, when we consider situations in which A and B are of different mag- 
nitude. 

Forest and Ruth (1990) give a fourth-order symplectic integrator with 6 substeps per step: 



S4B 



exp 



/tB\ 



- 



. 



exp 



\2c J 
/tA\ 



exp 



exp 



t) 

tB\ 



exp 



TB{l-k) 

2^ 



exp 



-TkA\ 
-^j.exp 



TB{l-k) 



2c 



\ c J 

= exp[rF + 0(r^)] 

where k = 2^'^ and c = 2 — k. Note that the middle 3 substeps move in the opposite direction to 
the integration as a whole. 

Higher-order integrators require progressively more substeps. Yoshida (1990) gives examples 
of 6th and 8th-order integrators using 14 and 30 substeps respectively. In the next section, we will 
show how to create what are effectively 4th and 6th order integrators (and in principle, 8th-order 
etc.) using fewer substeps than are required conventionally. 



3. Constructing Pseudo-Order Integrators 

Up to this point we have not considered the details of how H is split. Suppose that one part of 
the Hamiltonian is much smaller than the other, i.e. H = Ha + sHb, where e ^ 1. Now consider 
the error terms in the second-order integrator of equation (5): 



S2B 



exp 



rF+'-^[A,A,B] 



24 



[B,B,A] + --- 



One of the O^t"^) error terms is smaller than the other by a factor of e. 

Similarly, for the fourth order integrator: 

S4B = exp[rF + Oier^) + O(eV^) + 0{e^T^) + ©(eV^)] 

Some of these error terms are insignificant compared to others, and yet this was not taken into 
account when constructing the integrator. The only design criterion was that S4B should contain 
no error terms below the fifth power in the timestep. If we take into account the dependence of the 
error terms on both r and e, we can design more efficient symplectic integrators. 



To construct the new integrators, we again employ the BCH formula. Adapting the expression 
for the BCH formula given by Yoshida (1990), we have: 



exp(air^) • exp(6ieri3) 



exp 



{aiA + €biB)T + er^ 



aib 



aibf 



i) [A, B] + er^ (^) [A, A, B] + eV^ {^][B, B, A] 



+ e r ' 



H)[,,B,B,,,-.= (^)l.,AAABi-.v(^ 



[B,B,B,B,A] + 



where ai and 61 are constants. Additional fifth-order commutators are present; however, we will 
only require terms that contain either A o^ B once, since these are the type of error term we are 
seeking to eliminate. 

Applying the BCH formula twice, Yoshida (1990) gives an expression for a symmetric product 
of three exponential operators: 



exp(6ieri?) • exp(aiTA) • exp(6i eri?) 
exp 



{a,A + 2eh,B)T + er^ ( 4^^ [A, A, B] - eV^ (^ 



€T 



6 

^ ' ^ ^ [A,A,A, A, B] + e^T^ ' ^ 



[B,B,A] 



360 



360 



[B,B,B,B,A] + 



(6) 



Again we have neglected fifth-order terms that contain both A and B more than once. Note that 
there are no terms containing even powers of the timestep: Yoshida shows that this is a general 
property of any symmetric arrangement of exponential operators. From now on we will consider 
only symmetrical integrators because of this property. 

We need to extend equation (6) once more to get a pseudo-fourth order integrator, and twice 
more for a pseudo-sixth order one. By substituting 02^ for biB in equation (6), and substituting 
the righthand side of equation (6) for aiA, we get: 

exp(a2TA) • exp(6ieri?) • exp(air^) • exp(6ieri?) • exp(a2r^) 
= exp I (ai + 2a2)TA + IbierB + ^^^ (j) [("i + 2^2)^ - 602(01 + 02)] [A, A, B] 



+ eV3 ( ^A (4a2 - ai) [B, B, A] - er^ (^ ] [(ai + 202)^ - 30ai(ai + 02)^] [A, A, A, A, B] 



V36oy 



il6a2-7ai)[B,B,B,B,A] + --- 



(7) 



Finally, substituting the righthand side of equation (7) for aiA in equation (6), and replacing 
biB with 62 i3, we arrive at 



exp(62eri?) • exp(a2r^) • exp(6ieri?) • exp(airA) • exp(6i eri?) • exp(a2r^) • exp(62eT-B) 



6 



exp { (ai + 2a2)TA + 2(6i + b2)eTB + er^ 
ai + 2a2)ibi + b2f - 6a2bj 



{bi + 62) (ai + 202)^ - 60261(01 + 02) 



6 



[A,A,B] 



e\^ 



er 



6 



[B,B,A] 



jbi + 62)(qi + 202)"^ - 30a|6i(ai + 02) 
360 



[A,A,A,A,B] 



+ 6^5 



7(ai + 2a2)(bi + 62)^ - 60a2fcf (bi + 62)^ + 3002^1 
360 



[B,B,B,B,A] + 



(8) 



The first stage in converting these general expressions into specific integrators is to make the 
coefficients of the linear A and B terms equal to 1. This places two constraints on the values of 
the constants. We can then get what is effectively a 4th-order integrator by simply eliminating 
the [A, A, S] term from equation (7). The leading error terms will now be 0{e'^T'^) and 0{eT^). 
Provided that e is small enough, the largest error term will be 0{eT^), and the integrator effectively 
will be of fourth order in the timestep. Applying these conditions, we require 



ai + 2a2 = 1 

26i = 1 

1-602(1-02) = 

where we have used the first two of equations (9) in deriving the third. 



(9) 



Alternatively, we may construct an integrator in which each step begins by advancing Hb 
instead of -ff^. Unlike conventional symplectic integrators, such as 5*2^1 and S2B, we cannot 
use the same set of coefficients when exchanging A and B. Instead, we must derive a new set of 
coefficients by interchanging A and eB in equation (7) and then eliminating the new [A, A, B] term. 
When we do this, the first two of equations (9) remain as before, but the third expression becomes 



6a2 - 1 = 



(10) 



To get a pseudo-6th-order integrator, we eliminate terms containing [A, A, B] and [A, A, A, A, B]. 
This will produce an extra constraining equation, so we need an extra constant. We get this by 
using an integrator with the form of equation (8) instead of equation (7). The corresponding 
equations for the constants are 



oi + 2a2 = 


= 1 


2(61 + 62) = 


= 1 


1/2-60261(1-02) = 


= 


/2- 300^61(1-02)^ = 


= 



(11) 



If we prefer an integration step that begins by advancing Ha, we can interchange A and eB in 
equation (8), and eliminate the new [yl,yl, S] and [A, ^, ^, ^, i?] terms. In this case, the last two 



of equations (11) become 



1/4 - 6026? = 
7/16 - 15026? + 300261 = 



(12) 



The leading error terms for each of these integrators are 0{e^T^) and O^er'''). The latter will 
be dominant if e is small enough, so that the algorithms behave as 6th-order integrators. 



S4A* 



3.1. Pseudo-4th and 6th-Order Examples 

Solving equations (9) and (10), we obtain two pseudo-4th-order integrators: 



S4B* 



exp 



exp 



exp 



exp 



tA f 1 

^v "73 



exp y-^J . exp (^) . exp 



eTB\ 



2 V V3 



rF + eV ( ^-^^ [B,B,A] - ^[A,A,A,A,B] + 



24 ) ' ' ' ' 4320 

f2€TB\ (tA 



T) • ^"p \-y) • ^"p \r^) • ^^n 1" ' ■ '"p 



tA\ 

I • I 



eTB\ 



6^3 



er" 



rF + ^^[B,S,A] + ^^[AAAAi?] + 



where the asterisk in S4A* indicates that it only behaves as a 4th order integrator for certain values 
of r. 



Equations 11 and 12 give two pseudo-6th-order integrators: 



S6A* 



exp 



tA 



1 



^15 



• exp 



(heTB\ 



tA 



f5€TB\ 

it) 

3 



- 



• exp 



\2Vi5J 



/4eTB\ 



exp 



V 9 



• exp 






S6B* 



exp 

exp 
exp 
exp 



,,. + ,2^3/54-13^15^,, - ^1 .-^^^^^ 



648 






erB_ 

tA 



• exp 
1 



r^ 



1 

75 



5erS\ (tA\ ( heTB\ 



V5 



• exp 



^^2^3(13^5^5 

288 



eTB\ 

[S,S,A]+0(er^) 



Unlike the 4th-order algorithm of Forest and Ruth (1990) and the 6th-order integrators of 
Yoshida (1990), the algorithms above contain no substeps that move in the opposite direction to 



the main integration. An additional solution exists for each of equations (9), (11) and (12), however 
these have error terms with larger numerical coefficients than the integrators we show here. 

The same method can be used to generate a pseudo-8th-order integrator and so on. Each higher 
order will require only 2 more substeps than the previous one, since only one more commutator 
needs to be eliminated in each case. For example, to create a pseudo-8th-order integrator requires 
the elimination of the [A, A, A, A, A, A, B] term in addition to those that are absent from the 
pseudo-6th-order case. However, depending on the system to be integrated, there will come a point 
at which the e^r^ error term becomes the most important. In principle, one could devise another 
set of integrators that eliminates terms in e^r™" for small m, in addition to terms in er™. However, 
achieving each new order will generally require the elimination of more than one commutator term, 
so that these integrators increase in complexity much more rapidly than those described here. 

Murison and Chambers (1999) have independently derived the two 4th-order integrators above, 
among others, using a symbolic algebra package. Further results from that approach will follow 
in another paper. We note that the pseudo-order algorithms can be adapted to use indepen- 
dent timesteps for each planet (c./. (Saha and Tremaine 1994)), or to include close encounters 
((Duncan et al. 1998; Chambers 1999)). 



4. Numerical Comparisons 

In this section, we test the pseudo-4th and 6th order integrators derived in Section 3 against 
the well-known 2nd and 4th-order symplectic algorithms. We use the "mixed- variable" method 
of Wisdom and Holman (1991), in which the Hamiltonian is divided into a Keplerian part, Hk, 
and an interaction part. Hi. Under Hk^ each object moves on an unperturbed Keplerian orbit 
about the central body. Under i?/, each object remains fixed while receiving an impulse due to 
the gravitational perturbations of all the other objects except the central body. As suggested 
by Wisdom and Holman, we use Jacobi coordinates rather than barycentric coordinates. The 
integrations themselves were carried out using a modified version of the Mercury integrator package 
((Chambers and Migliorini 1997)). 

The pseudo-order integrators require that the ratio e = Hi/Hk ^1. In our first test, we 
integrate the orbits of the 4 inner planets of the solar system in the absence of the outer planets. 
In this case e ~ 10~^. Figure 1 shows the results of a 10000- year integration using the conventional 
2nd and 4th-order symplectic integrators, S2B and S4B, and the pseudo-order integrators S4B* 
and S6B*. For each integration, the maximum relative energy error is shown as a function of the 
step size. 

For the 2nd and 4th-order integrators, the maximum energy error is roughly proportional to 
r^ and r^ respectively, where r is the timestep. This is what we would expect to find. For the 
pseudo-4th and 6th-order integrators, the maximum energy error varies as r^ and r^. That is, they 
behave as 4th and 6th-order integrators, as we anticipated, despite the fact that they contain error 
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Fig. 1. — Maximum relative energy error versus step size for a 10000-year integration of the 4 
terrestrial planets using various symplectic integrators. 

Fig. 2. — Maximum relative energy error versus CPU time for a 10000-year integration of the 4 
terrestrial planets using various symplectic integrators. 
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Fig. 3. — Maximum relative energy error versus step size for a 10000-year integration of the 9 
planets using various symplectic integrators. 

Fig. 4. — Maximum relative energy error versus CPU time for a 10000-year integration of the 9 
planets using various symplectic integrators. 
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terms of lower order in the timestep. 

Using the mean relative energy error per integration instead of the maximum error gives results 
similar to Figure 1. The corresponding slopes are 2.10 it 0.05 for S2B, 3.9 it 0.3 for S4B, 4.6 it 0.3 
for S4B* and 6.4 ± 0.4 for S6B*. 

Figure 2 shows the amount of CPU time required for the integrations shown in Figure 1. For 
energy errors of 1 part in 10^ or 10^ there is not much to choose between the four algorithms. For 
higher levels of accuracy, S4B outperforms S2B. However, the pseudo integrators S4B* and S6B* 
do even better. At an accuracy of 1 part in lO^'^, they are roughly an order of magnitude faster 
than the conventional second-order integrator, and 3 times faster than the 4th-order integrator. 
For accuracies of better than 10^^^, S6B* shows greater performance than S4B*. 

The pseudo-4th order integrator is more efficient than the real 4th-order integrator for two 
reasons. It requires fewer substeps per time step, and it has a slightly smaller leading error term. 

As a more interesting test, we integrated the whole planetary system (Mercury to Pluto) for 
10000 years. Figure 3 shows the energy-error results of these integrations. The behaviour of S2B, 
S4B and S4B* is similar to that for the integrations of the terrestrial planets. However, the energy 
error for S6B* varies roughly as r^ rather than r^. It is not obvious why this should be, although 
the difference from the terrestrial-planet integration (Figure 1) is presumably due to the fact that 
e is two orders of magnitude larger in this case. 

Figure 4 shows the CPU time required for the integrations of the 9 planets. The results are 
similar to the integration of the inner planets, except that S6B* has only a marginal advantage 
over S4B* at the highest levels of accuracy. 

Since writing the original draft of this manuscript, we have become aware of the symplectic 
corrector method of Wisdom et al. (1996), which substantially improves the efficiency of the second- 
order symplectic integrator. We present the pseudo-order integrators as an alternative strategy for 
designing accurate algorithms. It is possible to devise other symplectic correctors using the same 
approach we use in Section 3 to design the integrator kernel: that is, by considering the dependence 
of the resulting error terms on e as well as r ((Mikkola 1997; Ranch and Holman 1999)). Finally 
we suggest that it may be possible to design symplectic correctors to improve the performance of 
pseudo-order algorithms, since the pseudo-order methods exhibit similar high-frequency oscillations 
in energy error to the second and 4th-order symplectic integrators (see Figure 5). 

In summary, we conclude that the new pseudo-order integrators outperform the widely-used 
2nd and 4th-order algorithms at all reasonable values of the energy error, for problems involving a 
dominant central mass. 

Research at Armagh Observatory is grant-aided by the Dept. of Education, Northern Ireland. 
The test integrations described in this paper were carried out using computers purchased on a 
PPARC research grant. 
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Fig. 5. — Relative energy error versus time for 10000-year integrations of the 9 planets using various 
symplectic integrators. 
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Figure Captions 

Figure 1: Maximum relative energy error versus step size for a 10000-year integration of tlie 4 
terrestrial planets using various symplectic integrators. 

Figure 2: Maximum relative energy error versus CPU time for a 10000-year integration of the 4 
terrestrial planets using various symplectic integrators. 

Figure 3: Maximum relative energy error versus step size for a 10000-year integration of the 9 
planets using various symplectic integrators. 

Figure 4: Maximum relative energy error versus CPU time for a 10000-year integration of the 9 
planets using various symplectic integrators. 

Figure 5: Relative energy error versus time for 10000-year integrations of the 9 planets using 
various symplectic integrators. 



